TK!2
=v
#1
:n

:s
b
:c
>>Sphere & cylinder interpenetration<<
#2
:n

:s
b
:c
--------------------------------------
#3
:n
R
:v
2.5
:s
i
:c
sphere radius
#4
:n
r
:v
1.5
:s
i
:c
cylinder radius
#5
:n
L
:v
3
:s
i
:c
distance of cyl axis & sphere center
#6
:n
b
:v
1.2472191289
:s
i
:c
upper limit of outer integral
#7
:n
Volume
:v
4.282590033853721
:s
o
:c
interpenetration volume
=u
=l
=f
#1
:n
dbintvl
:c
Double integration - composite Simpson's method
:t
Procedure
:i
fun,a,b,c,d,n
:o
VALUE
#2
:n
f
:t
Rule
:p
R
:a
x,y
:v
z
#3
:n
c
:t
Rule
:p
r,L
:a
x
:v
y
#4
:n
d
:t
Rule
:p
R
:a
x
:v
y
#1
:n
dbintvl
:z
" Notation:  fun    name of the rule or procedure function representing
"                   the integrand as a function of x,y
"            a,b    lower and upper limits of integration in x direction
"            c,d    names of the rule or procedure functions representing the"
"                   limits of integration in y direction as functions of x
"            n      number of integration intervals in x direction
"            VALUE  value of the double integral

" Definition:                               b   d(x)
" The procedure evaluates an integral      ô   ô
"                                          ³   ³  fun(x,y) dy dx 
"                                          õ   õ
"                                           a   c(x)
"
" using the composite Simpson's rule. The number of integration intervals
" in the x direction  n  must be given; the number of integration intervals
" in y direction  m  is evaluated so that m>=4 and hy<=hx. 

hx:= (b-a)/n
VALUE:= 0
x:= a
ki:= 4
for i=0 to n
   y:= apply(c,x)
   Hy:= apply(d,x) - y
   m:= 2*max(ceiling(Hy/hx/2),2)
   hy:= Hy/m
   value:= apply(fun,x,y) + apply(fun,x,apply(d,x))
   kj:= 2
   for j=1 to m-1
      y:= y+hy
      kj:= 6 - kj
      value:= value+kj*apply(fun,x,y) "accumulating the inner integral value
   next j
   value:= value*hy/3
   ki:= 6 - ki
   VALUE:= VALUE + ki*value           "accumulating the outer integral value
   if i=0 then VALUE:= VALUE - value
   x:= x+hx
next i
VALUE:= (VALUE-value)*hx/3
/E
E
#2
:n
f
:z
z^2 = abs(R^2 - x^2 - y^2)    " introduction of abs function prevents failure
                              " of sqrt when a small negative argument creeps
                              " in (for 0) due to floating point round-off
/E
E
#3
:n
c
:z
x^2 + (L-y)^2 = r^2     " equation of cylinder
/E
E
#4
:n
d
:z
x^2 + y^2 = R^2           " equation of sphere at z=0
/E
E
=r
#1
:r
"         >> D o u b l e   I n t e g r a t i o n  (variable limits) <<
:s
C
#3
:r
Volume = 4*dbintvl('f,0,b,'c,'d,10)
#4
:r
c(b) = d(b)
:s
C
#6
:r
" Notes ([PgDn] for more):
:s
C
#8
:r
" 1. Problem: Find the common volume occupied by a sphere with radius R=2.5
:s
C
#9
:r
"    and a cylinder with radius r=1.5 and axis at a distance L=3 from the
:s
C
#10
:r
"    center of the sphere. Suppose the sphere is centered at the origin of
:s
C
#11
:r
"    the coordinate system x,y,z, and the axis of the cylinder is parallel
:s
C
#12
:r
"    to the z-axis and lies in the  yz-plane. One quarter (by symmetry) of
:s
C
#13
:r
"    the interpenetration volume is given by
:s
C
#14
:r
"                                            b   d(x)
:s
C
#15
:r
"                                           ô   ô
:s
C
#16
:r
"                                           ³   ³  f(x,y) dy dx 
:s
C
#17
:r
"                                           õ   õ
:s
C
#18
:r
"                                            a   c(x)
:s
C
#19
:r
"
:s
C
#20
:r
"    where f(x,y) represents z evaluated from the equation of the sphere,
:s
C
#21
:r
"    and c(x) and d(x) represent respectively the y-coordinates (as func-
:s
C
#22
:r
"    tions of x) of the curves of intersection of the cylinder and the sphere
:s
C
#23
:r
"    with the xy-plane. The limits of the outer integral are a=0, and b is
:s
C
#24
:r
"    the x-coordinate of the intersection of y=c(x) and d(x).
:s
C
#25
:r
"
:s
C
#26
:r
" 2. The Function Procedure 'dbintvl' is used for evaluating the double
:s
C
#27
:r
"    integral by the composite Simpson's method. The functions f(x,y), c(x)
:s
C
#28
:r
"    and d(x) are defined in the associated rule functions 'f', 'c' and 'd'
:s
C
#29
:r
"    according to the description in 1.
:s
C
#30
:r
"
:s
C
#31
:r
" 3. There are two rules in this sheet. At the beginning, the first rule 
:s
C
#32
:r
"    is cancelled out and the Iterative Solver is used for evaluating the
:s
C
#33
:r
"    limit b (b=1 is a good guess). The resulting value of b is assigned as
:s
C
#34
:r
"    input, the second rule is cancelled out, and the first rule is used
:s
C
#35
:r
"    for evaluating the volume using the Direct Solver.
:s
C
#36
:r
"
:s
C
#37
:r
" 4. The Procedure Fuction 'dbintvl' may be simplified and the limit func-
:s
C
#38
:r
"    tions 'c' and 'd' may be omitted in cases of inner integrals with fixed
:s
C
#39
:r
"    limits (see the model DINTEGF). It also may be generalized for multiple
:s
C
#40
:r
"    integration.
:s
C
=p
=t
=n
=g
#1
:c
.000001
:i
10
:v
Yes
:l
60
:f
66
:m
0
:w
80
:b
Yes
:a
Yes
:r
Yes
:h
No
:t
1
:n
Yes
:d
LPT1
%Tr,f,0,1,0,1,0
